Giant Planet Engulfment by Evolved Giant Stars: Light Curves, Asteroseismology, and Survivability

Prior to contact between the star and planet, two effects bring the planet and the stellar surface closer together: stellar radius expansion and orbital decay. Orbital decay is driven by tidal friction, most likely due to turbulent dissipation in the star's convective envelope (e.g., Zahn 1977, 1989; Vick & Lai 2020). Drag forces exerted by the stellar corona and wind are negligible (Duncan & Lissauer 1998).

When the planet makes contact with the star, drag forces begin to dominate over tidal friction. The "grazing" hydrodynamical interaction of the star and planet is complex and three-dimensional (see MacLeod & Loeb 2020a; Yarza et al. 2022). Various observable phenomena, such as the expulsion of stellar matter (MacLeod & Loeb 2020b; Lau et al. 2022) and shock-powered optical and X-ray transients (Metzger et al. 2012; Stephan et al. 2020), may occur during the grazing phase. Those are beyond the scope of this study. We focus on the later "inspiral" phase of engulfment, when the planet is completely immersed in the envelope.

We use a set of fiducial stellar models obtained with MESA-r22.05.1 (see Section 5 for software details), meant as representative models of RGB and AGB stars. We evolved a nonrotating star of initial mass 1.50 M_⊙ from the pre-MS through the end of the thermally pulsing AGB stage. We adopted an initial chemical composition using the protosolar abundances of Asplund et al. (2009). We included wind-driven mass loss using the prescriptions of Reimers (1975) from the zero-age MS (ZAMS) through the RGB and Blöcker (1995) on the AGB, with respective scaling factors η_R = 0.5 and η_B = 0.1. We treated convective boundaries using the Schwarzschild criterion and neglected overshooting. Our stellar models are snapshots of this single evolutionary sequence. Table 1 summarizes their major properties. Each model is labeled in the form ABCxyz, where ABC refers to the star's evolutionary stage and xyz is its approximate radius in units of R_⊙.

Table 1. Properties of the Fiducial Host Star Models

Model	R⋆	M⋆	${M}_{\mathrm{core}}$	Egrav	Etot	L⋆	Teff	∣λ∣	τdyn	τKH	${ \mathcal N }(3{M}_{{\rm{J}}})$
	(R⊙)	(M⊙)	(M⊙)	(erg)	(erg)	(L⊙)	(K)	⋯	(yr)	(yr)	⋯
RGB50	50.1	1.493	0.362	−2.1 × 1047	−6.9 × 1046	549	3950	1.8	0.015	2500	150
RGB100	100.0	1.420	0.424	−1.1 × 1047	−2.4 × 1046	1511	3599	2.2	0.042	420	340
RGB150	149.3	1.357	0.474	−7.1 × 1046	−9.9 × 1045	2672	3396	3.1	0.079	150	550
AGB200	203.7	1.278	0.555	−4.7 × 1046	−3.5 × 1045	4124	3242	4.9	0.13	60	833
AGB275	276.7	1.000	0.568	−1.8 × 1046	+1.5 × 1045	5602	3003	4.0	0.23	20	1950

Note. All stellar models were evolved from a 1.50 M_⊙ ZAMS model. Definitions: ${M}_{\mathrm{core}}$ is the enclosed mass at the bottom of the convection zone. The global dynamical time is ${\tau }_{\mathrm{dyn}}\equiv {({R}_{\star }^{3}/{{GM}}_{\star })}^{1/2}$ and Kelvin–Helmholtz time is ${\tau }_{\mathrm{KH}}\equiv {{GM}}_{\star }^{2}/({R}_{\star }{L}_{\star })$. The gravitational binding energy of the envelope is E_grav. The total binding energy E_tot and factor λ are defined in Equation (21). The total number of orbits undergone by a 3 M_J giant planet is ${ \mathcal N }(3{M}_{{\rm{J}}})$.

Download table as:  ASCIITypeset image

The drag force acting on an engulfed planet is determined by the flow of stellar matter around it. We assume that the planet instantaneously follows a circular orbit with radius a and velocity ${\boldsymbol{v}}={v}_{{\rm{K}}}(a){\hat{{\boldsymbol{e}}}}_{\phi }$, with

$$\begin{eqnarray}&&{v}_{{\rm{K}}}={\left(\displaystyle \frac{{{GM}}_{r}}{r}\right)}^{1/2},\end{eqnarray} \tag{ 1 }$$

where M_r is the stellar mass enclosed by a shell of radius r. The planet's Mach number is

$$\begin{eqnarray}&&{ \mathcal M }={v}_{{\rm{K}}}/{c}_{s},\end{eqnarray} \tag{ 2 }$$

where c_s is the local speed of sound in the stellar envelope. We can safely assume that v_K and c_s are much greater than the bulk velocity of fluid near the planet due to stellar rotation and convection. The planet's motion is always supersonic, with $1\lesssim { \mathcal M }\lesssim 10$ in general and $1\lesssim { \mathcal M }\lesssim 2$ in the deepest portion of the envelope (where the most energy is dissipated). The run of ${ \mathcal M }$ for stellar model AGB200 is shown by the solid curve in Figure 1.

Zoom In Zoom Out Reset image size

Despite the supersonic motion through the envelope, the interior of the inspiraling giant planet remains in hydrostatic equilibrium for most of the inspiral because the ram pressure exerted by the medium is small compared to the planet's central pressure. The drag on the planet comprises hydrodynamic friction (a.k.a., ram pressure or geometric friction), which arises from the nonuniform fluid stress exerted across the object's surface, and gravitational dynamical friction, which arises from the wake formed behind the object by gravitational focusing. The importance of each component is determined by the planet's radius R_p and gravitational focusing (a.k.a., Bondi–Hoyle) radius

$$\begin{eqnarray}&&{R}_{{\rm{a}}}=\displaystyle \frac{2{{GM}}_{{\rm{p}}}}{{v}^{2}},\end{eqnarray} \tag{ 3 }$$

where M_p is the planet's mass. For an object with velocity v moving through a uniform medium of density ρ, the hydrodynamic and gravitational drag forces may be written

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{\rm{d}}}=-{C}_{{\rm{d}}}\pi {R}_{{\rm{p}}}^{2}\rho v{\boldsymbol{v}}\propto {v}^{2},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{\rm{g}}}=-{C}_{{\rm{g}}}\pi {R}_{{\rm{a}}}^{2}\rho v{\boldsymbol{v}}\propto \displaystyle \frac{1}{{v}^{2}}.\end{eqnarray} \tag{ 5 }$$

The coefficients C_d and C_g (see below) are determined by the global flow structure, which depend on the dimensionless quantities ${ \mathcal M }$ and R_a/R_p. Note that our definition of C_d absorbs a leading factor (1/2) that others prefer to be explicit. All else being equal, gravitational (hydrodynamical) drag dominates when v is small (large), which tends to occur at large (small) r within the envelope. The transition between the drag regimes occurs when v is of the order of the planet's surface escape velocity, ${v}_{{\rm{e}}}={(2{{GM}}_{{\rm{p}}}/{R}_{{\rm{p}}})}^{1/2}$.

Thun et al. (2016) and Yarza et al. (2022) studied the flow around a gravitating supersonic projectile using 2D axisymmetric and 3D hydrodynamics simulations, respectively. Both used a "wind tunnel" setup to characterize the steady-state flow structure and numerically derive the hydrodynamic and gravitational drag forces on the projectile (see also MacLeod et al. 2017). For supersonic motion, a shock develops ahead of the projectile. The structure of the postshock material depends on the Mach number and the "compactness" parameter R_a/R_p. Thun et al. (2016) found that the stand-off distance R_sh of the shock from the center of the projectile is given approximately by

$$\begin{eqnarray}&&{R}_{\mathrm{sh}}={R}_{{\rm{p}}}\max \left(1,\eta \right),\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&\eta \equiv \displaystyle \frac{g(\gamma )}{2}\displaystyle \frac{{{ \mathcal M }}^{2}}{{{ \mathcal M }}^{2}-1}\displaystyle \frac{{R}_{{\rm{a}}}}{{R}_{{\rm{p}}}},\end{eqnarray} \tag{ 7 }$$

is a nonlinearity parameter (Kim & Kim 2009) and the factor g(γ) depends on the adiabatic exponent γ. We take γ = 5/3 and g(5/3) = 1 throughout this work. The hydrodynamic crossing time (∼R_sh/v) is short compared to the Keplerian orbital period and the orbital decay time, so we assume that the shock structure is in a steady state determined by local conditions. Figure 1 shows the run of R_a/R_p and η for model AGB 200 for two planet masses, M_p = 1 M_J and 10 M_J.

When gravitational focusing is negligible (R_a ≪ R_p), the postshock material is compressed ahead of the projectile and rarefied behind it. The total drag force is then dominated by hydrodynamic friction. The drag coefficient is a weak function of Mach number and can be estimated using laboratory ballistics data (Bailey & Hiatt 1972). We use the following fitting formula

$$\begin{eqnarray}&&{C}_{{\rm{d}}}=0.375+0.125\tanh [1.75({ \mathcal M }-1)],\end{eqnarray} \tag{ 8 }$$

which provides an adequate estimate for a nongravitating projectile with Reynolds number ≫1 (appropriate for our case).

In the limit of strong gravitational focusing (R_a ≫ R_p), postshock material accumulates around the projectile, creating a quasi-spherical "halo" of subsonic matter with radius ≈ R_a (Thun et al. 2016). This structure arises from the impenetrable boundary condition at the projectile's surface; it does not occur in simulations with absorbing or outflowing boundary conditions (e.g., MacLeod et al. 2017; Li et al. 2020). The halo is approximately in hydrostatic equilibrium and exerts almost uniform pressure on the projectile; thus hydrodynamic drag is negligible (C_d ≃ 0) in this limit. The gravitational drag coefficient is given by Ostriker (1999) and Thun et al. (2016)

$$\begin{eqnarray}&&{C}_{{\rm{g}}}=\mathrm{ln}\left(\displaystyle \frac{{R}_{\max }}{{R}_{\mathrm{sh}}}\right)-\displaystyle \frac{1}{2}\mathrm{ln}\left(\displaystyle \frac{{{ \mathcal M }}^{2}}{{{ \mathcal M }}^{2}-1}\right),\end{eqnarray} \tag{ 9 }$$

where ${R}_{\max }$ is the maximum extent of the wake formed behind the sphere.

In order for the "wind tunnel" results of Thun et al. (2016) to be applicable to the motion of a body orbiting within a stratified stellar envelope, ${R}_{\max }$ must not exceed the local density scale height, H_ρ . The importance of the density gradient is measured by the parameter

$$\begin{eqnarray}&&{\epsilon }_{\rho }=\displaystyle \frac{{R}_{\mathrm{sh}}}{{H}_{\rho }}.\end{eqnarray} \tag{ 10 }$$

Yarza et al. (2022) conducted wind tunnel simulations with 0.1 ≤ _ρ ≤ 1 and 0.3 ≤ R_p/R_a ≤ 1. They found that the presence of a density gradient introduces asymmetry to the flow structure around the object but otherwise does not change the main conclusions of Thun et al. (2016). Notably, the total drag force increases by less than a factor of ∼2 relative to that exerted by a uniform medium. For our cases, we find that _ρ ≲ 0.2 for r ≲ 0.95 R_⋆, indicating that corrections to the local drag force law due to the envelope's density gradient are negligible for giant planets whose "spheres of influence" (∼few R_sh) are well separated from the stellar surface.

The results summarized above indicate that, because the structure of the flow around the planet depends mainly on R_a/R_p, the hydrodynamic and gravitational regimes of the drag force are almost mutually exclusive. We adopt the following prescription for the total drag force during the inspiral phase

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{drag}}=-{\hat{{\boldsymbol{e}}}}_{\phi }\max \left[{C}_{{\rm{d}}}\pi {R}_{{\rm{p}}}^{2},{C}_{{\rm{g}}}\pi {R}_{{\rm{a}}}^{2}\right]\rho {v}_{{\rm{K}}}^{2},\end{eqnarray} \tag{ 11 }$$

where C_d is given by Equation (8) and C_g by Equation (9) with ${R}_{\max }={H}_{\rho }$.

The rate of orbital decay is directly related to the rate of nonconservative work by drag forces

$$\begin{eqnarray}&&\displaystyle \frac{\dot{a}}{a}=\displaystyle \frac{{{\boldsymbol{F}}}_{\mathrm{drag}}\cdot {\boldsymbol{v}}}{| {E}_{\mathrm{orb}}| }\equiv -\displaystyle \frac{1}{{\tau }_{\mathrm{insp}}},\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&{E}_{\mathrm{orb}}=-\displaystyle \frac{{{GM}}_{r}(a){M}_{{\rm{p}}}}{2a},\end{eqnarray} \tag{ 13 }$$

and τ_insp is a characteristic inspiral time at a given separation. We assume that the orbital energy loss is deposited in the vicinity of the planet as heat. The heat deposition rate is

$$\begin{eqnarray}&&{L}_{\mathrm{drag}}=-{{\boldsymbol{F}}}_{\mathrm{drag}}\cdot {\boldsymbol{v}}=-{\dot{E}}_{\mathrm{orb}}.\end{eqnarray} \tag{ 14 }$$

Heat is deposited at the shock that develops around the planet (Thun et al. 2016; Yarza et al. 2022) and diffuses via convection into the surroundings (e.g., Wilson & Nordhaus 2022). For further discussion of the role of stellar convection, see Section 4.

When gravitational dynamical friction dominates, the inspiral timescale is roughly

$$\begin{eqnarray}&&{\tau }_{\mathrm{insp}}\sim \displaystyle \frac{1}{{C}_{{\rm{g}}}}\displaystyle \frac{{M}_{r}}{{M}_{{\rm{p}}}}{\tau }_{\mathrm{orb}},\end{eqnarray} \tag{ 15 }$$

where ${\tau }_{\mathrm{orb}}=2\pi {({a}^{3}/{{GM}}_{r})}^{1/3}$ is the local Keplerian orbital period. When hydrodynamic drag dominates, it is roughly

$$\begin{eqnarray}&&{\tau }_{\mathrm{insp}}\sim \displaystyle \frac{1}{{C}_{{\rm{d}}}}\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{r}}{\left(\displaystyle \frac{a}{{R}_{{\rm{p}}}}\right)}^{2}{\tau }_{\mathrm{orb}}.\end{eqnarray} \tag{ 16 }$$

In all cases we consider, the planet undergoes ${ \mathcal N }\sim {10}^{2}$–10³ orbits in total during the inspiral.

3.2.1. Angular Momentum Deposition

Drag forces also exert a torque that reduces the planet's orbital angular momentum and spins up the star. In principle, one could study the differential rotation induced by deposition of the planet's angular momentum in successive layers of the star using MESA. Since our focus is the effect of energy deposition by the planet, we neglect spin-up in our calculations. It suffices to say that the engulfment would greatly increase the star's rotation rate (e.g., Livio & Soker 2002; Carlberg et al. 2009, 2012; Privitera et al. 2016; Stephan et al. 2020). Assuming the star spins up rigidly and has negligible initial angular momentum, the postengulfment rotation rate is

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\star }=\displaystyle \frac{1}{{k}_{\star }}\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{\star }}{{\rm{\Omega }}}_{{\rm{K}}\star },\end{eqnarray} \tag{ 17 }$$

where Ω_K⋆ is the Keplerian angular velocity at the stellar surface (a.k.a., the breakup rate) and where the stellar moment of inertia is ${k}_{\star }{M}_{\star }{R}_{\star }^{2}$ with k_⋆ ≈ 0.1. For M_p = 1–10 M_J and M_⋆ ≈ M_⊙, final rotation rates ≈ 0.01–0.1 Ω_K⋆ are expected. Such rapid rotation would be a lasting signature of engulfment detectable via spectroscopic broadening, photometric modulation, and asteroseismic mode splitting. It could also generate magnetic dynamo activity in the envelope (e.g., Livio & Soker 2002; Nordhaus & Blackman 2006).

Engulfed planets are usually destroyed in the stellar interior. This can occur either catastrophically or gradually, depending on the planet's properties and the envelope conditions.

Livio & Soker (1984) proposed that an engulfed substellar body is thermally disrupted ("dissolved") upon reaching a location where the local sound speed exceeds the body's surface escape velocity. Though frequently quoted in the literature (e.g., Siess & Livio 1999a, 1999b; Carlberg et al. 2009; Aguilera-Gómez et al. 2016; Privitera et al. 2016; Cabezón et al. 2022; Lau et al. 2022), this criterion is only the first step needed for dissolution. The second is to compare the inspiral time with the timescale for heating of the planetary interior by the ambient medium. For irradiated giant planets, the rate-limiting processes are conductive and radiative heat transfer, as the rising entropy of the outer layers suppresses convection (Guillot et al. 1996; Arras & Bildsten 2006). This implies a heating time far exceeding the predicted inspiral time; hence, we do not consider thermal disruption to be relevant for giant planets.

Alternatively, planets can be disrupted by ram pressure or ablation due to their supersonic motion in a stellar envelope (Jia & Spruit 2018). Neither of these is relevant for us because the planet is much denser than its surroundings, even near the base of the convection zone.

Based on these considerations, we assume that an engulfed planet survives in the stellar envelope until it fills its Roche lobe, whereupon it is undergoes rapid tidal disruption (see Reyes-Ruiz & López 1999; Nordhaus & Blackman 2006; Nordhaus et al. 2011; Guidarelli et al. 2022). The orbital radius a_dis where this occurs is (Eggleton 1983)

$$\begin{eqnarray}\begin{array}{rcl}{a}_{\mathrm{dis}} & \simeq & 2{R}_{{\rm{p}}}{\left(\displaystyle \frac{{M}_{\mathrm{core}}}{{M}_{{\rm{p}}}}\right)}^{1/3}\\ & \approx & 1.8{R}_{\odot }\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{{\rm{J}}}}\right){\left(\displaystyle \frac{{M}_{\mathrm{core}}}{0.5{M}_{\odot }}\displaystyle \frac{{M}_{{\rm{J}}}}{{M}_{{\rm{p}}}}\right)}^{1/3},\end{array}\end{eqnarray} \tag{ 18 }$$

where ${M}_{\mathrm{core}}$ is the mass of the compact stellar core. Planets are disrupted long before making physical contact with the core (r ≲ 0.1 R_⊙). We do not consider the subsequent evolution of tidal debris in this work. However, previous studies have considered various possibilities (Siess & Livio 1999a, 1999b; Reyes-Ruiz & López 1999; Nordhaus & Blackman 2006; Nordhaus et al. 2011; Guidarelli et al. 2022).

When a planet undergoes inspiral from an initial separation a₁ ≃ 1 R_⋆ to a final separation a₂ ≪ a₁, the total energy deposited in the envelope is given by

$$\begin{eqnarray}&&W=\left|-\displaystyle \frac{{{GM}}_{r}({a}_{2}){M}_{{\rm{p}}}}{2{a}_{2}}+\displaystyle \frac{{{GM}}_{\star }{M}_{{\rm{p}}}}{2{a}_{1}}\right|\simeq \displaystyle \frac{{{GM}}_{\mathrm{core}}{M}_{{\rm{p}}}}{2{a}_{2}},\end{eqnarray} \tag{ 19 }$$

setting the overall energy budget for the inspiral. If heat deposition ceases (or at least slows down) when the planet is tidally disrupted, then by setting a₂ = a_dis, we find an estimate of the characteristic energy budget

$$\begin{eqnarray}\begin{array}{rcl}{W}_{\mathrm{insp}} & = & \displaystyle \frac{G{({M}_{\mathrm{core}}^{2}{M}_{{\rm{p}}}^{4})}^{1/3}}{4{R}_{{\rm{p}}}}\\ & \approx & 5.6\times {10}^{44}\,\mathrm{erg}\\ & & \times \,{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{{\rm{J}}}}\right)}^{-1}{\left(\displaystyle \frac{{M}_{\mathrm{core}}}{0.5{M}_{\odot }}\right)}^{2/3}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{{\rm{J}}}}\right)}^{4/3}.\end{array}\end{eqnarray} \tag{ 20 }$$

This should be compared with the total binding energy of the stellar envelope

$$\begin{eqnarray}&&{E}_{\mathrm{tot}}={\int }_{{M}_{\mathrm{core}}}^{{M}_{\star }}\left(-\displaystyle \frac{{{GM}}_{r}}{r}+u\right){\rm{d}}{M}_{r}\equiv -\displaystyle \frac{{{GM}}_{\star }{M}_{\mathrm{env}}}{\lambda {R}_{\star }},\end{eqnarray} \tag{ 21 }$$

where ${M}_{\mathrm{env}}={M}_{\star }-{M}_{\mathrm{core}}$ is the mass of the envelope and u is the specific internal energy of the stellar gas. Although the disruption radius a_dis varies as a function of planet mass, the mass enclosed by the terminal orbit is always very close to ${M}_{\mathrm{core}}$. Table 1 lists the values of E_tot and λ for our fiducial stellar models. We also give the gravitational binding energy E_grav of the envelope, computed by omitting u in Equation (21).

Both ∣E_tot∣ and ∣E_grav∣ range from 10⁴⁵ to 10⁴⁷ erg in order of magnitude. The envelope's E_tot is usually negative (λ > 0) but becomes positive (λ < 0) on the upper AGB, where ionization energy overwhelms gravitational potential energy (Paczyński & Ziółkowski 1968).

Comparison of W_insp and ∣E_tot∣ gives a sense of whether the deposition of heat by the planet leads to significant changes in the overall structure of the stellar envelope during inspiral. By setting W_insp = ∣E_tot∣ and solving for M_p, we obtain a characteristic planetary mass for which we expect a major structural adjustment

$$\begin{eqnarray}\begin{array}{rcl}{M}_{{\rm{p}},\mathrm{crit}} & = & {\left(\displaystyle \frac{4}{| \lambda | }\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\star }}\right)}^{3/4}\displaystyle \frac{{({M}_{\star }{M}_{\mathrm{env}})}^{3/4}}{{M}_{\mathrm{core}}^{1/2}}\\ & \approx & 2.8{M}_{{\rm{J}}}{\left(\displaystyle \frac{4}{| \lambda | }\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{{\rm{J}}}}\displaystyle \frac{200{R}_{\odot }}{{R}_{\star }}\right)}^{3/4}\\ & & \times \,{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\displaystyle \frac{{M}_{\mathrm{env}}}{0.5{M}_{\odot }}\right)}^{3/4}{\left(\displaystyle \frac{0.5{M}_{\odot }}{{M}_{\mathrm{core}}}\right)}^{1/2}.\end{array}\end{eqnarray} \tag{ 22 }$$

However, as we show in Section 5, energetic arguments alone cannot predict the range of possible outcomes.

Our calculation of the inspiral energy budget is somewhat incomplete in that we ignore the evolution of the debris from the planet's tidal disruption (see previous section). The continued release of gravitational potential energy by the sinking debris (e.g., Siess & Livio 1999a, 1999b) would contribute additional energy W_sink ∼ W_insp over the timescale τ_sink, assuming that the debris sinks to a final radius ∼ r_tide/2. Since the central temperature of a gas giant or brown dwarf is typically much less than the ambient temperature, ∼T(r_tide), the debris would also absorb energy (−W_abs) as heat from the environment over the timescale τ_abs (Harpaz & Soker 1994). The quantity of heat absorbed is

$$\begin{eqnarray}\begin{array}{rcl}| {W}_{\mathrm{abs}}| & \sim & \displaystyle \frac{{M}_{{\rm{p}}}}{{m}_{p}}{k}_{{\rm{B}}}T({r}_{\mathrm{tide}})\\ & \approx & 2\times {10}^{44}\,\mathrm{erg}\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{{\rm{J}}}}\right)\left(\displaystyle \frac{T({r}_{\mathrm{tide}})}{{10}^{6}\,{\rm{K}}}\right),\end{array}\end{eqnarray} \tag{ 23 }$$

where m_p is the mass of a proton and k_B is Boltzmann's constant. This happens also to be of the order of W_insp in the cases we consider. These additional process act as residual sources and sinks of energy after a planetary inspiral. The timescales τ_sink and τ_abs are somewhat uncertain; they depend in detail on the dynamical evolution of the tidal debris. If either timescale is comparable to, or shorter than, the inspiral time just before tidal disruption (τ_late; see Section 3.6), these effects could potentially alter the evolution of the host star to a significant degree. On the other hand, if τ_abs ≈ τ_sink, then the effects may cancel one another because W_insp ∼ −W_orb. It would be possible to incorporate these effects in MESA using a parameterized prescription for the sinking and thermal evolution of tidal debris. We leave this for future work.

A crude but useful approximation of the inspiral trajectory can be obtained by assuming that an engulfed planet has no effect on the stellar envelope. We call this the "passive envelope" approximation. As the inspiral is brief compared to the stellar evolutionary timescale, a snapshot of the stellar structure is sufficient to compute the planet's trajectory.

The right-hand side of Equation (12) determines the inspiral rate $\dot{a}$ in terms of the properties of the envelope and planet. In the passive envelope approximation, the planet's orbital separation a at time t is given by

$$\begin{eqnarray}&&t-{t}_{1}=-{\int }_{a}^{{a}_{1}}\displaystyle \frac{{\rm{d}}r}{\dot{a}(r;{M}_{{\rm{p}}},{R}_{{\rm{p}}})},\end{eqnarray} \tag{ 24 }$$

where a₁ is the separation of the planet at time t₁. The instantaneous heating rate due to drag (L_drag) along the planet's trajectory can also be computed.

As the planet spirals inward, the cumulative heat deposited in the star grows like 1/a. Roughly half of the total heat is deposited between r = 2 a_dis and r = a_dis. We dub this the "late" stage of inspiral. This stage is critical for determining the qualitative behavior of the stellar envelope response. We denote the energy deposited in the envelope during this stage as ${W}_{\mathrm{late}}\simeq {W}_{\max }/2$. The duration of the late inspiral, denoted τ_late, is of the order of the local inspiral time τ_insp evaluated at r = 2 a_dis; in the passive envelope approximation, τ_late can be computed exactly by way of Equation (24). The average heating rate due to drag during the late inspiral is L_late ≡ W_late/τ_late.

In Figure 2, we show L_late as a function of M_p for all stellar models, as computed under the passive envelope approximation. For a given star, L_late often exceeds the intrinsic luminosity L_⋆ and is nearly proportional to M_p. The latter fact implies that τ_late depends weakly on M_p, in qualitative agreement with the scaling ${\tau }_{\mathrm{insp}}\propto {M}_{{\rm{p}}}{a}_{\mathrm{dis}}^{2}\propto {M}_{{\rm{p}}}^{1/3}$ predicted by Equations (16) and (18). Our MESA experiments confirm that the passive envelope approximation predicts the planet's trajectory accurately in many cases; they also allow for a prediction of when the approximation breaks down.

Zoom In Zoom Out Reset image size

The MESA equation of state (EOS) is a blend of the OPAL (Rogers & Nayfonov 2002), SCVH (Saumon et al. 1995), FreeEOS, ⁷ HELM (Timmes & Swesty 2000), PC (Potekhin & Chabrier 2010), and Skye (Jermyn et al. 2021) EOSs. Radiative opacities are primarily from OPAL (Iglesias & Rogers 1993, 1996), with low-temperature data from Ferguson et al. (2005). Electron conduction opacities are from Cassisi et al. (2007) and Blouin et al. (2020). Nuclear reaction rates are from JINA REACLIB (Cyburt et al. 2010), NACRE (Angulo et al. 1999), and additional tabulated weak reaction rates (Fuller et al. 1985; Oda et al. 1994; Langanke & Martínez-Pinedo 2000). Screening is included via the prescription of Chugunov et al. (2007). Thermal neutrino loss rates are from Itoh et al. (1996). The inlist and extension files required to reproduce our results are available online. ⁸

In all runs, we enable MESA's 1D hydrodynamics module and treated convection using the recently added TDC option (Jermyn et al. 2023). We include localized heating from an engulfed planet using the drag force formula from Section 3 and the shellular approximation described in Section 4. We evolve the the planetary orbit simultaneously with the stellar model using Equation (12) with the initial a = 0.95 R_⋆. We use Equation (18) to determine when the planet is tidally disrupted and switch off the extra heating immediately thereafter. In order to focus on the hydrodynamical effects of inspiral heating, we do not allow for wind-driven mass loss during these simulations.

By default, we terminated each MESA run after an elapsed time ≈10 τ_KH since the planet's disruption, where τ_KH is the stellar Kelvin–Helmholtz time. In most cases, this is ample time for the star to revert to its original structure. However, a few runs ended early because convergence could not be attained with a reasonable simulation time step. In those cases, planetary engulfment had triggered supersonic expansion of the envelope (see Section 5.4).

We refer to our MESA experiments using an abbreviated notation: for example, RGB150-3MJ refers to the simulation with stellar model RGB150 (see Table 1) and a 3 M_J planet.

Figures 4(a), (b), and (c) show the simultaneous evolution of the planet and stellar structure with stellar models RGB50, RGB100, and RGB150 for M_p = 1–10 M_J. Figure 4(d) shows AGB200-1MJ, 3MJ, and 5MJ, and Figure 5 shows AGB275-1MJ and 3MJ. Each shows the evolution of the planetary orbital radius, stellar radius, and stellar bolometric luminosity and absolute magnitude. To facilitate comparison between simulations, we define the time at which the planet is tidally disrupted to be t = 0, with negative and positive values of t corresponding to the time before and after disruption.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The inspiral of a massive planet triggers expansion and brightening of the envelope, followed by protracted Kelvin–Helmholtz contraction and gradual dimming. The more massive the planet, the larger the disturbance it creates. For M_p = 3–10 M_J, the star's initial brightening is plausibly detectable with a ground-based optical or near-infrared telescope, with typical photometric deviations of >0.1 mag over ≈1 yr. The effective temperature T_eff also decreases by ∼10–500 K as the star expands, so a mild amount of intrinsic reddening is expected. We emphasize that the brightening is powered by the planet's descent in the deep interior, not by a "grazing" interaction at the stellar surface. A star may not brighten noticeably for decades to centuries after engulfment begins.

The expected shape of the light curve in this phase depends mostly on the host star's evolutionary stage: compact, early-stage red giants such as RGB50 and RGB100 tend to brighten in a monotonic fashion before leveling off. The puffy, more-evolved RGB150, AGB200, and AGB275 display two distinct peaks in their light curves separated by ≈1 yr. Smaller secondary peaks can also be seen for RGB50 and RGB100. We discuss the physical origin of the double-peaked light curve below.

The results for AGB275 are notable because the star's convective envelope has E_tot > 0, perhaps making it more susceptible to ejection (Paczyński & Ziółkowski 1968). Although planets of 1 M_J and 3 M_J cause relatively large disturbances, they do not unbind the envelope, likely because of efficient radiative cooling in the expanding outer layers (see below).

We now examine a subset of our MESA experiments in which planetary engulfment disrupts the stellar envelope, namely AGB200-10MJ, AGB275-5MJ, and AGB275-10MJ. These represent a transitional regime between the quasi-static and subsonic responses described above and the dynamical ejection of matter during a true common envelope event.

Figure 6 shows the results of experiment AGB275-5MJ, including the evolution of T_eff. The evolution of the stellar radius resembles a more extreme version of experiment AGB275-3MJ in that the star undergoes large-scale expansion and contraction, reaching a maximum radius of 1500 R_⊙ (≈7 au). However, this representation obscures the rich hydrodynamics of interior mass shells (see below). The light curve and T_eff evolution reflect this complexity more effectively. The star initially brightens by ≈1.5 mag over ≈2 yr before fading by ≈4 mag on a similar timescale. At the same time, it reddens from T_eff ≈ 3000 K to just under 1000 K. When the expansion reverses, the star abruptly returns almost to the same peak luminosity and T_eff ≈ 2000 K. Unlike in the previous experiments, the star's outer layers are essentially in free fall at this stage. The luminosity at t ≳ 4 yr is generated by an accretion shock at the interface between the expanding/collapsing "ejecta" and the quasi-static interior, rather than by residual heat from the planetary inspiral.

Zoom In Zoom Out Reset image size

The final "spike" in the light curve, accompanied by an increase in T_eff to ≈5500 K, coincides with the reaccretion of the photosphere and is analogous to a shock breakout in a stellar explosion. Its duration is d/v_ff, where v_ff is the free fall velocity at the shock and d is the thickness of an outer shell of optical depth ${\bar{\tau }}_{r}=c/{v}_{\mathrm{ff}}\sim {10}^{4}$.

Despite the drastic expansion of the envelope, no material is ejected from the star above the escape velocity. We note that MESA's EOS takes the ionization energy of hydrogen and helium into account. This extra energy reservoir is evidently insufficient to eject the envelope in our scenario. However, we have not accounted for the formation of dust grains as the envelope expands and cools. Dust-driven winds may play a crucial role in the late stages of a common envelope event (Glanz & Perets 2018). If the grain opacity is large enough during the envelope's initial expansion, it may reduce radiative losses and allow the outer layers to reach escape velocity. The evolution predicted by MESA at late times (t ≳ 5 yr) therefore should be regarded with caution.

Other effects not included in our simulations, such as stellar rotation and magnetic activity, would likely influence the evolution of a disrupted, dusty stellar envelope. If the spinning-up of the envelope by the planet generates a magnetic dynamo, then magnetic "cool spots" on the stellar surface may form dust grains more readily; this would likely change the geometry and mass-loss rate of a dust-driven outflow (e.g., Soker 1998b; Livio & Soker 2002; Nordhaus & Blackman 2006; Rapoport et al. 2021).

Figure 7 shows AGB200-10MJ and AGB275-10MJ. In these cases, the MESA runs were terminated early because the required time step for convergence became prohibitively short. Consequently, we cannot characterize the ejecta dynamics or light curve at late times in these cases. However, each model was evolved through the planet's tidal disruption with satisfactory accuracy. Broadly speaking, the star's initial expansion phase resembles that of AGB275-5MJ in both cases.

Zoom In Zoom Out Reset image size

5.4.1. Hydrodynamics of the Interior

Figures 6 and 7 show that the stellar photosphere eventually stops expanding and falls back. In Figure 8, we show selected snapshots of the radial velocity and entropy profiles of the stellar interior from AGB275-5MJ. The envelope does not expand and contract monotonically throughout the evolution. Instead, a given mass shell undergoes multiple phases of expansion and contraction in general, and different mass shells can expand and contract simultaneously. The maximum expansion velocity is ≈40% of the escape speed.

Zoom In Zoom Out Reset image size

Initially, the star expands and contracts uniformly in the sense that the sign of the radial velocity is the same everywhere. The expansion is initially adiabatic, but subsequently the entropy of the outer layers drops abruptly as a result of hydrogen recombination (see Section 5.4.2). Radiative energy transport temporarily dominates in a large fraction of the envelope, due to the presence of an inverted entropy gradient (Figure 8). Radiative cooling in these regions further deprives the expanded envelope of pressure support, leading to collapse. This accounts for the precipitous decline in stellar luminosity after the initial broad peak. However, the stellar interior eventually becomes overpressurized from the collapse and rebounds. Thus, for a time, the star contains two dynamically disconnected regions: a subsonic expanding interior and an envelope of ejecta in supersonic free fall. A strong shock is present at the contact layer, visible as a discontinuity in the velocity and entropy profiles. This shock converts the kinetic energy of collapse into thermal and ionization energy, allowing the star to return to equilibrium.

Although the planet fails to eject the stellar envelope in our MESA models, we have nonetheless identified a significant transition in their qualitative behavior: for very advanced stages of stellar evolution, there is a critical planetary mass above which the outer layers of the envelope expand supersonically and develop shocks (even if they do not become unbound). For AGB200 (AGB275), the critical mass is between 5 and 10 M_J (3 and 5 M_J). This agrees roughly with the quantity M_p,crit estimated in Section 3.4 via energy budget considerations (Equation (22)).

However, the inspiral energy budget alone cannot determine whether the stellar model undergoes a nonlinear hydrodynamical response. This depends rather on the rate of heat deposition during the planet's late inspiral (L_late, see Section 3.6), which is proportional to M_p for a given star. Assuming that the star initially undergoes quasi-static evolution, with L_late giving the work per unit time against self-gravity, the expansion rate at the surface, ${v}_{\exp }$, can be estimated as follows

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{late}} & \sim & \displaystyle \frac{{{GM}}_{\star }{M}_{\mathrm{env}}}{{R}_{\star }^{2}}{v}_{\exp },\\ \Longrightarrow {v}_{\exp } & \sim & \displaystyle \frac{{L}_{\mathrm{late}}{R}_{\star }^{2}}{{{GM}}_{\star }{M}_{\mathrm{env}}}\\ & \approx & 6\,\mathrm{km}\,{{\rm{s}}}^{-1}\left(\displaystyle \frac{{L}_{\mathrm{late}}}{{10}^{5}{L}_{\odot }}\right){\left(\displaystyle \frac{{R}_{\star }}{200{R}_{\odot }}\right)}^{2}\\ & & \times \,{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\displaystyle \frac{{M}_{\mathrm{env}}}{0.5{M}_{\odot }}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 30 }$$

This can be compared with the sound speed near the surface

$$\begin{eqnarray}&&{c}_{s}={\left(\displaystyle \frac{\gamma {k}_{{\rm\small{B}}}T}{{\mu }_{\mathrm{gas}}}\right)}^{1/2}\approx 12\,\mathrm{km}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{3\gamma }{5}\displaystyle \frac{T}{{10}^{4}\,{\rm{K}}}\displaystyle \frac{{m}_{p}}{{\mu }_{\mathrm{gas}}}\right)}^{1/2},\end{eqnarray} \tag{ 31 }$$

where γ and μ_gas are the gas' adiabatic index and average mass per particle. In AGB stars, ${v}_{\exp }$ can be comparable to c_s because R_⋆ is large. If the expansion continues with L_late held constant, ${v}_{\exp }$ increases in proportion to ${R}_{\star }^{2}$ and c_s drops due to adiabatic cooling. Consequently, the expansion becomes supersonic for large-enough L_late (∝M_p), as we have seen.

We have verified this picture by conducting additional AGB200-10MJ runs in which the drag force on the planet is artificially reduced by a constant factor. Although the same amount of heat is deposited in the star in each case, the stellar response becomes quasi-static when the drag force is reduced by roughly an order of magnitude. This is because (a) the envelope expanded more slowly and (b) convection carries a greater fraction of L_late.

5.4.2. Recombination and Radiative Losses (without Dust)

The interior structure profiles reveal that a large portion of the envelope undergoes recombination as a result of adiabatic cooling during the star's initial expansion. The top panel of Figure 9 shows the evolution of μ_gas throughout the envelope for experiment AGB275-5MJ as a proxy for the predominant phase of hydrogen. As the envelope expands following the planet's inspiral, more than half of the envelope's mass undergoes recombination. (The temperature near the surface can be low enough to form molecules and dust grains, but this is a secondary effect in terms of energetics.)

Zoom In Zoom Out Reset image size

In the study of common envelope evolution, the role of ionization energy in ejecting the envelope is uncertain (see the discussions in the reviews of Ivanova et al. 2013 and Röpke & De Marco 2023, as well as, e.g., Sabach et al. 2017; Grichener et al. 2018; Ivanova 2018). The crux of the issue is whether a sufficient portion of the ionization energy can be converted to mechanical work to accelerate the envelope to escape velocity. Several previous studies, including Sabach et al. (2017), Grichener et al. (2018), and Wilson & Nordhaus (2019), have argued that convection can efficiently transport heat derived from orbital and ionization energy from the deep interior to radiative zones near the surface, reducing the amount of energy available to do work (but see Ivanova 2018). We find that much of the envelope's initial ionization energy is lost as radiation during the initial large-scale expansion of the envelope (see Section 5.4.2). Subsequently, the remaining heat deposited by the planet flows outward, reionizing the envelope in the process. This implies that the fluid absorbs a substantial amount of thermal energy that could otherwise have been used as work to eject the envelope. The ionization energy therefore acts as a "buffer" for the planet's orbital energy, arguably preventing envelope ejection rather than aiding it. These findings demonstrate the nontrivial effects of ionization energy transport in common envelope evolution: even if it does not contribute to ejecting the envelope, it may be necessary to include ionization effects in numerical calculations in order to predict a system's evolution and observable characteristics accurately. However, we reiterate that our conclusions are provisional, due to the lack of dust effects in our simulations.

Recombination has two main effects on stellar matter in our MESA models: the ionization energy converts to radiation and the opacity plummets. Even though the envelope remains optically thick, the reduced opacity enables efficient radiative cooling. To demonstrate this, the middle and bottom panels of Figure 9 show the evolution of the radiative opacity κ reported by MESA and the local radiative cooling timescale τ_rad. The cooling time is defined as

$$\begin{eqnarray}&&{\tau }_{\mathrm{rad}}={\bar{\tau }}_{r}\displaystyle \frac{{R}_{\star }-r}{c}\displaystyle \frac{{P}_{\mathrm{gas}}}{{P}_{\mathrm{rad}}},\end{eqnarray} \tag{ 32 }$$

where ${\bar{\tau }}_{r}={\int }_{r}^{{R}_{\star }}\kappa (r^{\prime} )\rho (r^{\prime} )\,{\rm{d}}r^{\prime}$ is the optical depth at radius r, c is the speed of light, and P_gas and P_rad are the local gas pressure and radiation pressure, respectively. In the layers that undergo recombination, τ_rad becomes comparable to, or shorter than, the envelope expansion time ${\tau }_{\exp }\approx 1\,\mathrm{yr}$. Roughly speaking, layers with ${\tau }_{\mathrm{rad}}\lt {\tau }_{\exp }$ can cool efficiently.

This line of reasoning also reproduces the maximum luminosity during the initial expansion of the envelope. Assuming that a shell of mass M_sh loses all of its hydrogen ionization energy as radiation over a time ${\tau }_{\exp }$, we find that the star's maximum luminosity is given by

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\max } & = & \displaystyle \frac{{{XM}}_{\mathrm{sh}}{q}_{{\rm{H}}}}{{\tau }_{\exp }}\\ & \approx & 3.2\times {10}^{4}{L}_{\odot }\left(\displaystyle \frac{X}{0.75}\displaystyle \frac{{M}_{\mathrm{sh}}}{0.2{M}_{\odot }}\right){\left(\displaystyle \frac{{\tau }_{\exp }}{1\,\mathrm{yr}}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 33 }$$

where X is the hydrogen mass fraction and q_H = 13.6 eV/m_p . This agrees with the MESA results for AGB200-10MJ, AGB275-5MJ, and AGB275-10MJ. (All else being equal, accounting for the ionization energy of helium increases ${L}_{\max }$ by 10%.) Recombination is also responsible for the first peak in the light curves of experiments RGB150-10MJ, AGB200-5MJ, and AGB275-3MJ (Figures 4(c) and (d) and 5, respectively). In those cases, the envelope expands by a smaller factor and remains subsonic, so M_sh and ${L}_{\max }$ are reduced relative to the reference values above.

As noted above, dust grains are expected to form as the stellar envelope expands and cools. Their additional opacity may reduce radiative energy losses, allowing more efficient ejection of the envelope. In order for this to occur, the grain opacity must be large enough that τ_rad remains longer than ${\tau }_{\exp }$. Referring to the bottom panel of Figure 9, we see that the grain opacity must be 2–3 orders of magnitude greater than the low-temperature MESA opacity to eliminate radiative losses completely in AGB275-5MJ.

Previous studies have suggested that planet engulfment (or tidal disruption) may lead to an observable optical/infrared transient, including Retter & Marom (2003), Retter et al. (2006, 2007), Bear et al. (2011a), Metzger et al. (2012), Kashi & Soker (2017), Metzger et al. (2017), MacLeod et al. (2018), Soker (2018), Kashi et al. (2019), Stephan et al. (2020), and Gurevich et al. (2022). Our results confirm this prediction and characterize the relationship between transient properties and underlying planetary and stellar parameters.

The most informative quantity to constrain the engulfed planet's mass is the excess stellar luminosity. Let ΔL(t) be the difference between the star's instantaneous bolometric luminosity L_tot(t) and its intrinsic luminosity L_⋆. Let the maximum value be ${\rm{\Delta }}{L}_{\max }$ occurring at time ${t}_{\max }$. The upper panel of Figure 10 shows the quantity ${\rm{\Delta }}{L}_{\max }/{L}_{\star }$ for our MESA simulations. For stars that experience quasi-static or subsonic responses (which generally correspond to ${\rm{\Delta }}{L}_{\max }/{L}_{\star }\lesssim 1$), it may be possible to estimate the planetary mass from the peak luminosity, provided that the host's properties are well constrained prior to the flare-up. However, when the planet creates a major disturbance in the star (Section 5.4), nonlinear effects make the relation between M_p and ${\rm{\Delta }}{L}_{\max }$ more complex.

Zoom In Zoom Out Reset image size

The maximum change of the host's intrinsic color (which we characterize via T_eff) also correlates with planetary mass. The effective temperature tends to decrease somewhat during the transient, reddening the star. The lower panel of Figure 10 shows the maximum change of T_eff from its initial value. Cases with more-massive planets and more-evolved host stars exhibit larger T_eff changes. Again, however, nonlinear effects dominate when the planet disrupts the envelope.

The brightening and fading timescales of engulfment-powered transients are also important. In Figure 11, we show the values of the empirical quantities Δt_bri and Δt_fade, defined as ${\rm{\Delta }}{t}_{\mathrm{bri}}={t}_{\max }-{t}_{{\rm{A}}}$ and ${\rm{\Delta }}{t}_{\mathrm{fade}}={t}_{{\rm{B}}}-{t}_{\max }$, where ${\rm{\Delta }}L({t}_{{\rm{A}}})={\rm{\Delta }}L({t}_{{\rm{B}}})=0.1{\rm{\Delta }}{L}_{\max }$ and ${t}_{{\rm{A}}}\lt {t}_{\max }\lt {t}_{{\rm{B}}}$.

Zoom In Zoom Out Reset image size

For systems with quasi-static or subsonic responses, the fading time matches the Kelvin–Helmholtz time

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{KH}} & = & \displaystyle \frac{{{GM}}_{\star }^{2}}{{R}_{\star }{L}_{\star }}\\ & \approx & 40\,\mathrm{yr}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{2}{\left(\displaystyle \frac{{R}_{\star }}{200{R}_{\odot }}\displaystyle \frac{{L}_{\star }}{4000{L}_{\odot }}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 34 }$$

This shows that the star contracts quasi-statically as the planet's orbital energy is lost as radiation. Since τ_KH is longer than a century for stars with radii ≲100 R_⊙, the secular dimming of a star following engulfment is measurable only for tip of the red giant branch (TRGB) and AGB stars. However, it may be possible to detect the contraction of stars in the ≈100–150 R_⊙ range by indirect means (see below).

The brightening time varies within the range of ∼1–5 yr across almost all of our runs, without any apparent trend with planetary mass for a given stellar model or between stellar models. It is long compared to the global dynamical timescale

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{dyn}} & = & {\left(\displaystyle \frac{{R}_{\star }^{3}}{{{GM}}_{\star }}\right)}^{1/2}\\ & \approx & 0.14\,\mathrm{yr}{\left(\displaystyle \frac{{R}_{\star }}{200{R}_{\odot }}\right)}^{3/2}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-1/2},\end{array}\end{eqnarray} \tag{ 35 }$$

which reflects the star's expansion below the escape velocity. If we assume that the star undergoes quasi-static expansion during the planet's late inspiral, we estimate the expansion timescale (see Equation (30))

$$\begin{eqnarray}&&{\tau }_{\exp }=\displaystyle \frac{{{GM}}_{\star }{M}_{\mathrm{env}}}{{R}_{\star }{L}_{\mathrm{late}}}\sim {\tau }_{\mathrm{KH}}\left(\displaystyle \frac{{L}_{\star }}{{L}_{\mathrm{late}}}\right).\end{eqnarray} \tag{ 36 }$$

Using L_late ≈ 10⁴–10⁶ L_⊙ for the appropriate stellar parameters (Figure 2), Equation (36) reproduces the range of Δt_bri seen in our MESA results. Run RGB50-1MJ is an exception because the stellar response is so weak that small acoustic oscillations (with period ≈ τ_dyn) excited during the late inspiral dominate the light curve.

6.1.1. Weak Red Transients

Broadly speaking, the electromagnetic signatures of planetary inspiral can be separated into two qualitative categories based on whether the stellar response was quasi-static/subsonic or supersonic. The main signature of a quasi-static/subsonic envelope is an abrupt (≈1 yr) increase in luminosity with mild reddening (see Figures 4 and 5 for examples). TRGB and AGB stars may exhibit a prominent "double peak" in bolometric luminosity due to H recombination. In most cases, the transient phase is followed by a long plateau (10–100 yr) at nearly constant L and T_eff as Kelvin–Helmholtz contraction takes over. The bolometric amplitudes of these transients are ∼0.1–1 mag for planets between 3 and 10 M_J. Despite their modest amplitudes, these "weak red transients" may be observable by ground- and space-based wide-field time-domain surveys operating at optical and near-infrared wavelengths, such as the Zwicky Transient Facility, the Vera C. Rubin Observatory, and the Nancy Grace Roman Space Telescope.

We estimate the rate of weak red transients in the Galaxy based on the rate of WD formation (Γ_WD ∼ 1 yr⁻¹) and the occurrence fraction of warm Jupiters around FGK dwarfs (f_WJ ≈ 0.03; Cumming et al. 2008; Mayor et al. 2011)

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{WRT}}={f}_{\mathrm{WJ}}{{\rm{\Gamma }}}_{\mathrm{WD}}\approx 0.03\,{\mathrm{yr}}^{-1}\left(\displaystyle \frac{{f}_{\mathrm{WJ}}}{0.03}\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{WD}}}{1\,{\mathrm{yr}}^{-1}}\right)\end{eqnarray} \tag{ 37 }$$

.

6.1.2. Red Novae from Engulfment on the AGB

Despite some physical uncertainties and computational limitations, our MESA experiments clearly predict that sufficiently massive planets can disrupt an AGB star's envelope, producing a major eruption lasting several years (Figures 6 and 7). These eruptions likely resemble the luminous red novae (LRNe) produced by coalescing binary stars (e.g., Tylenda & Soker 2006). However, they are much dimmer at peak brightness and evolve more slowly than typical LRNe (e.g., Kochanek et al. 2014; Karambelkar et al. 2023). Additionally, their progenitors are already bright, dusty, late-type sources by virtue of their evolutionary stage, whereas the progenitors of LRNe display a wide range of properties.

If the fraction of WD progenitors with a cold Jupiter is f_CJ, the rate of engulfment-powered RNe in the Galaxy is

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{RNe}}\approx 0.1\,{\mathrm{yr}}^{-1}\left(\displaystyle \frac{{f}_{\mathrm{CJ}}}{0.1}\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{WD}}}{1\,{\mathrm{yr}}^{-1}}\right),\end{eqnarray} \tag{ 38 }$$

where we again estimate f_CJ based on the occurrence rate among Sun-like stars (Fernandes et al. 2019; Fulton et al. 2021). If the giant planet occurrence rate is greater for stars somewhat more massive than the Sun (Johnson et al. 2010; Reffert et al. 2015; Jones et al. 2016; Ghezzi et al. 2018), then a somewhat larger value of f_CJ may be appropriate for typical single WD progenitors (M_⋆ ≈ 1.5–3 M_⊙ on the MS). The higher rate and peak luminosity of these eruptions relative to weak red transients makes them a more-promising target population for transient surveys.

When an AGB star is disrupted by an engulfed planet, it reaches a typical peak luminosity M_bol ≈ −6 and lasts several years. These properties, combined with the expected event rate ∼0.1 yr⁻¹, suggest that planetary engulfment events could be responsible for a significant fraction of the lower-luminosity RNe observed in the Galaxy (Kochanek et al. 2014; Howitt et al. 2020). The transient OGLE-2002-BLG-360 (Tylenda et al. 2013), which was "less violent" and of longer duration than a typical LRN, and whose progenitor was a late-type giant star, may be an example.

6.1.3. Observational Complications

In the "weak red transient" regime, the dimming timescale following planet engulfment is longer than the brightening timescale by a factor of ∼10–10³, depending on evolutionary stage. Thus, there are ∼10–10³ times more stars in the dimming phase than the brightening phase at any given time. For the most part, the dimming is too gradual to be measured on human timescales. However, it may be possible to detect the star's quasi-static contraction using asteroseismology. Specifically, a contracting star would display a secular increasing trend in the large seismic frequency spacing ${\rm{\Delta }}\nu \propto 1/{\tau }_{\mathrm{dyn}}\propto {R}_{\star }^{-3/2}$. The Kepler mission produced a homogeneous set of Δν measurements with a typical fractional uncertainty of 10⁻³ for several thousand RGB stars monitored continuously over ≈4 yr (Yu et al. 2018). Stars contracting at rates $| {\dot{R}}_{\star }| /{R}_{\star }\gtrsim {10}^{-3}\,{\mathrm{yr}}^{-1}$ may have a detectable Δν trend.

Additionally, giant stars that have engulfed a giant planet rotate more rapidly than is typically observed or expected based on theoretical models of single star evolution (see Section 3.2.1 and references therein). Rapid rotation would cause rotational mode splitting in the seismic spectrum and would be a longer-lasting signature of planet engulfment than brightness changes. The spin-down timescale of a red giant due to stellar mass loss is $\simeq 0.1{({M}_{\star }/{\dot{M}}_{{\rm{w}}})({R}_{\star }/{r}_{{\rm{w}}})}^{2}$, where ${\dot{M}}_{{\rm{w}}}$ is the rate of mass loss due to a stellar wind and r_w is the radius from which the wind is launched. For r_w ≈ R_⋆, the spin-down time is ≈10% of the star's remaining lifetime ($\sim {M}_{\star }/{\dot{M}}_{{\rm{w}}}$). A giant displaying both strong mode splitting and an upward secular drift of Δν would be a strong candidate for having recently engulfed a substellar companion.

For completeness, we note that an engulfed planet or brown dwarf can excite stellar oscillations directly during its inspiral phase (e.g., Soker 1992a, 1992b; Gagnier & Pejcha 2023). This has possible ramifications for wind-driven mass loss and the detectability of ongoing engulfment events.

As mentioned previously, this work is motivated in part by the question of whether an engulfed giant planet can survive on a short-period orbit around a WD after ejecting the stellar envelope (Vanderburg et al. 2020; Chamandy et al. 2021; Lagos et al. 2021; Merlov et al. 2021). In all of our MESA experiments, the planet undergoes Roche-lobe overflow (RLO) before the end of the simulation. We have assumed that, once RLO begins, the planet is disrupted rapidly compared to the ongoing dynamical evolution of the envelope. This assumption bears further scrutiny.

In principle, RLO of a short-period giant planet can lead to stable mass transfer, halting or even reversing its orbital decay (Valsecchi et al. 2014; Jackson et al. 2016; Jia & Spruit 2017). In our case, the inspiral/mass-transfer timescale is much shorter than the planet's thermal timescale. The planet's radius therefore evolves at a fixed entropy per unit mass, increasing or remaining nearly constant (e.g., Zapolsky & Salpeter 1969; Paxton et al. 2013). Thus, the mass transfer is unstable and runaway disruption ensues.

RLO can be prevented if the stellar envelope expands to such a degree that the planet's inspiral time becomes longer than the global dynamical time, "stalling" the planet at a > a_dis until the envelope has dissipated completely. In "successful" common envelope ejections, these self-regulating behaviors prevent a complete merger and set the post-common envelope binary separation (e.g., Ivanova et al. 2013; Clayton et al. 2017; Gagnier & Pejcha 2023). The only experiment in which the inspiral became somewhat self-regulating was AGB275-10MJ (Figure 7). However, this merely delayed RLO by ≈10 yr, much less than the expected timescale of dust-driven mass loss from an AGB or post-common envelope (e.g., Glanz & Perets 2018). This suggests that planetary survival is possible only for M_p ≳ 10 M_J and only when engulfment occurs during the thermally pulsing AGB stage. This, in turn, implies that short-period giant planets can only be found orbiting WDs with carbon/oxygen cores (as opposed to close "WD+brown dwarf" binaries, which often contain helium-core WDs; see Zorotovic & Schreiber 2022 and references therein).

Currently, two intact, short-period giant planets (or planet candidates) are known orbiting single WDs. The transiting planet WD 1856+534 b (Vanderburg et al. 2020) orbits a (0.58 ± 0.04) M_⊙ WD at a separation of 0.02 au (≈4 R_⊙). Its mass is constrained to be 0.84 M_J < M_p < 13.8 M_J (Vanderburg et al. 2020; Xu et al. 2021). Our MESA simulations suggest that a planet would struggle to avoid tidal disruption through most of the allowed mass range, even if it succeeded in ejecting the envelope with the aid of dust-driven winds. For this system, high-eccentricity tidal migration is a plausible alternative formation scenario (Muñoz & Petrovich 2020; O'Connor et al. 2021; Stephan et al. 2021). Meanwhile, WD 0141–675 has an astrometric planet candidate with mass ${9.3}_{-1.1}^{+2.5}{M}_{{\rm{J}}}$ and semimajor axis 0.17 au (≈34 R_⊙) according to the Gaia DR3 astrometric solution (Gaia Collaboration et al. 2022). If confirmed, this system would be difficult to explain as a common envelope survivor: the planet's large orbit implies an insufficient energy budget to disrupt the envelope (but see Bear & Soker 2011, 2014; Bear et al. 2021; Chamandy et al. 2021). The microlensing giant planet MOA-2010-BLG-477L b is also associated with a WD (Blackman et al. 2021). Based on its large sky-projected separation from the host (2.8 ± 0.5 au), this planet likely avoided engulfment during late-stage stellar evolution.

Over the years, several studies have reported short-period planet candidates orbiting horizontal branch stars or hot subdwarfs (e.g., Geier et al. 2009; Setiawan et al. 2010; Charpinet et al. 2011). However, none has been confirmed to date (e.g., Norris et al. 2011; Jones & Jenkins 2014; Krzesinski 2015). Such objects would have survived engulfment during the first RGB and may have ejected a portion of the stellar envelope (e.g., Bear & Soker 2011, 2014). We do not find conditions under which this can occur for M_p ≤ 10 M_J. Some recent works have suggested that planetary survival is possible if engulfment coincides with the host's core helium flash (Bear et al. 2011b, 2021; Merlov et al. 2021). Our MESA extension could be used to reexamine this scenario.